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] Abstract. A general semi-analytical method for accurate and efficient numerical 

f^ calculation of the dielectrically screened ( "dressed" ) potential around a non-relativistic 

*7^ test particle moving in an isotropic, collisionlcss, unmagnetised plasma is presented. 

£h The method requires no approximations and is illustrated using results calculated for 

c/3 two cases taken from the MSc thesis of the first author: test particles with velocities 

i ^ above and below the ion sound speed in plasmas with Maxwellian ions and warm 

Qh electrons. The idea that the fluctuation spectrum of a plasma can be described as a 

c/2 superposition of the fields around non-interacting dressed test particles is an expression 

of the quasiparticle concept, which has also been expressed in the development of the 



oscillation-centre and pseudo-orbit formalisms. 
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1. Introduction 

The Fokker-Planck-Landau equation incorporating dielectric screening now known as 
the Balescu-Lenard equation was also derived by Thompson and Hubbard [H El E] from 
simpler statistical physics arguments. In the Thompson-Hubbard approach the diffusion 
coefficient was calculated from a fluctuation spectrum obtained by superimposing the 
dielectrically screened fields of independently moving particles, an approach called the 
dressed test particle picture by Rostoker jl] . 

In this picture, unperturbed "test particles" replace the actual particles in the 
plasma, though in reality the trajectories are perturbed slightly by the fluctuations, 
giving rise to the linear dielectric response and velocity-space diffusion. The dressed 
test particle picture was developed in Fourier, (u>, k), representation rather than in 
real space-time, (x, t), and thus it was difficult to visualise the actual nature of the 
screened potential surrounding each particle. This created interest in inverting the 
Fourier transform to better understand Balescu-Lenard kinetic theory. 

The first chapter of the MSc thesis of the first author [5] was on relativistic plasma 
response functions [DJ. However Chapters 2 and 3 of the thesis were restricted to the 
non-relativistic case and calculated the screened potential of a dressed test particle in 
real space, using both asymptotic approximation methods and "exact" calculations in 
which the triple Fourier integral was evaluated in a way that exploited a highly efficient 
numerical algorithm for calculating a special function r](z) defined for the purpose. 

Much of the two relevant chapters of reference [5] has recently been published more- 
or-less verbatim [7], but details of the numerical approach were omitted. Our main 
purpose in this paper is to present the mathematical analysis and numerical method 
used in reference [5] in the hope they may prove useful in teaching and further research. 

A survey [5] of the '60s literature on test-particle screening is reprinted in reference 
[?]. The problem was revisited occasionally in the '70s, e.g. [8], and '80s, but in the 
last decade and a half there has been a considerable revival of interest in this topic in 
the context of dusty plasmas. Some of these more recent papers are mentioned briefly 
below. 

Ishihara and Vladimirov [DJ applied screened potentials to charged dust particles 
in a plasma. They calculated the potential in the wake of a moving test particle and 
showed it contained periodic minima. This can result in an attractive force between 
dust particles, providing a mechanism for the formation of Coulomb crystals. Other 
studies have considered the effect of Landau damping [TD1 fTTj . neutral collisions [ID] . 
magnetic fields [12] and the wake of a dipole [13]. Lapenta [H] considered a derivation 
of the screened potential in real space, rather than Fourier space, with the aim of 
including nonlinear effects in the wake. Numerical simulations have played an important 
role in validating analytic results and finding potentials beyond the point particle 
approximation for objects such as rods [15]) and multiple particles [16j [17]. 

While particle-in-cell simulations [THJ [TB] allow more physics to be included, 
including nonlinearity [19] [20], the approaches based on linear dielectric response, 
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e.g. [HI El HU [16], in principle allow more resolution and accuracy (within the linear 
regime) and should be numerically less intensive. However, to simplify the problem the 
dielectric response function is usually approximated and/or integrations are performed 
only approximately. 

An exception appears to be the work of Lampe and Joyce [TDJ [TB], who used 
Fast Fourier Transforms to perform the Fourier inversions without the need to make 
approximations. However, because the inversions were performed only once and 
tabulated to avoid the need for recalculation, calculational details and timings were 
not given in these papers. 

In this paper we present the efficient numerical approach [5] that allowed an 
extensive numerical study even using the limited computer power available in the mid- 
'60s. The method does not require the dielectric response function to be approximated, 
so it includes Landau damping accurately, and was applied to dielectric constants 
corresponding to both Lorentzian and Maxwellian plasma distribution functions. We 
have found that the method allows a full, accurate calculation of the wake structure to 
be computed in a few minutes on a modern laptop. 

Details of the numerical techniques used and the special function rj were omitted 
from reference [7j so they are, until now, unpublished. In section 15] we review the 
Fourier representation of dielectric screening and present the strategy introduced in 
reference [5] for performing the 3-dimensional integral over k in spherical polars, by 
evaluating the infinite integral over k first. This is the most difficult part of the 
numerical Fourier inversion as numerical integration of rapidly oscillatory integrands is 
notoriously difficult, but this difficulty is circumvented by drawing on the mathematical 
tools of complex analysis and the literature on special functions to develop an efficient 
algorithm for calculating this integral in terms of a special function, rj(z), thus effectively 
performing the fc-integral analytically. 

In section [3] we introduce the coordinate system used to perform the remaining 2- 
dimensional integral over solid angle numerically, and in section [4] we indicate how this is 
implemented in a reconstruction of the code, using gf ortran [2T]. Some of the previous 
cases studied [3 [7] are recalculated and replotted as validation of the reconstructed 



code and a few new, but related plots are presented as well. Appendix A gives details 
of the special function r]{z) and its relation to the exponential integral E\{z) and the 
auxiliary function for sine and cosine integrals, f(z). A subroutine for calculating rj(z) 
is provided online as supplementary material. 

In section [5] we attempt to trace how the patterns of thought developed from the 
MSc studies of the first author have influenced the course of his career. 
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2. Dressed test particles 

A test particle with charge q moving through a plasma at constant, non-relativistic 
velocity v produces the following dielectrically screened potential [7] (in SI units) 

q f d 3 k expfzk • (x — v £) — A|k|] 

ip = 7 x™oJ J2irf Pe(k-v ,k) ' () 

where Sq is the permittivity of free space and e(u, k) is the frequency, u), and wavevector, 
k, -dependent plasma dielectric constant. In an isotropic, collisionless, unmagnetised 
plasma the dielectric constant is of the form 

e( W ,k) = e(w,|k|) = l + ^^, (2) 

where §(u/k) (called the polarisation function in references [5j U\) will be discussed more 
explicitly below. For now all we need to assume is that it obey the reality condition, 
$(-z/) = $*(z/), and the stability condition [7j, Re[$(z/)] 1/2 > 0, for all real finite 
z/, where [^(z/)] 1 / 2 = |$(z/)| 1//2 exp |iarg$(z/). In dlh , A is a regularisation parameter 
required to interpret the integral as k = |k| — > oo. 

Writing r = x — vo£, we have the potential in the rest frame of the test particle 

tfr) = S- lim / ^fc ^(«k-r-A|k|) 
y ' e \->+aJ {2ixf fc 2 + $(k-v ) 

where k is the unit vector in the direction of k. A key insight in reference j5] was that 
the 3-dimensional integral in pi) could be reduced to the 2-dimensional integral below 
by transforming to spherical polars and integrating over k using the special function 



r){-) defined in (A.l ). Making use of the identity (A. 3) differentiated twice with respect 



to a = k-r + iX, with (3 = [$(k • v )] 1/2 = a/$, we find 

¥ ,( r ) = _ i l. i im / ^ffivfT/YCk.r + zA)^ , (4) 

where dQ(k) is an element of solid angle such that d 3 k = k 2 dQ(k) and S 2 is the unit 



sphere. Using the identity (A. 13), which gives rj"(z) = r](z) — 1/z, we can also write <p 



in the form of the "bare" Coulomb potential (po(r), 

Mr) " 4^ ' (5) 

plus a correction term, which may be interpreted as the potential of the screening charge 
"dressing" the particle, 

y>(r) = Mr) -i±L a ^ v^7/(k.rv^) . (6) 

The first term on the RHS of ([6]), the bare potential fo(r), comes from the 1/z 
term in r]"(z), which is divergent as z — > 0. We used the Plemelj formula to interpret 
lim A ^ +0 l/(k-r + zA) as P/k-r — 27r<5(k-r), where V denotes the principal part operator 
and 5(-) is the Dirac delta function. (The integral over the principal part term vanishes, 
so only the delta-function term contributes to the potential.) We have set A = in 
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the second term, the dressing potential, as, from (A. 8), the singularity in rj( k-rv$) at 
k-r = is sufficiently weak that regularisation is not required. 

The function $(•), as derived in the standard way from the linearised Vlasov 
equation, is given by 

• (£)=E<£ /"*-#-■ m 

\kj V p J-™ u/k-v 

Here co ps denotes the plasma frequency, {e 2 n s /eQm s ) 1 ^ 2 , for species s, with n s the 
unperturbed number density, m s the mass, e s the charge, and g s (v) the one-dimensional 
projection of the velocity distribution function / s (v). For a non-relativistic Maxwellian 
plasma, 

/ ^2 V 



$J^x)=k 2 



h ' '' ^ 



V2xF ■'' ' 



1- V2xF I -^= I +iy-xexp 



(8) 



i-Ds 

where /cds = (e^ns/eo^s) 1 '' 2 is the inverse Debye length for species s, T s being the 
temperature in energy units, and F(x) the Dawson function, defined by [22] 

F{Q=e-? f C e t2 dt. (9) 



The RHS of (I8J can also be written as — (k 2 :)s /2)Z'(x/\/2) where Z(Q is the Plasma 
Dispersion Function 



2.1. Limitations of the method 

It should be noted that the method requires that /c 2 e(k-vo, k) be of the form k 2 + ^(k) 



in order for the identity (A. 3) to be used to evaluate the infinite integral over k. This 
limits the applicability of the method to classical, collisionless, unmagnetised plasmas 
as explained below. The specific form we have assumed for $ in (J2J) and ([7]) is further 
limited to isotropic plasmas, but the method would also apply to the anisotropic case, 
complicating only the integration over solid angle. 

The classical polarisation function $ involves numerators, k-(df s /dv), and 
denominators, u — k-v, arising from solving the linearised Vlasov equation, while in 
the quantum case [21] the numerators are replaced by {m s /%) [/ s (v + Kk/2m s )) — / s (v — 
hk/2m s )], which reduces to the classical form as % — > 0. To obtain the form in ([7]), 
or its anisotropic generalisation, we divide both numerators and denominators by k, so 
that all ^-dependences of the numerators are removed and the numerators depend on k 
only through the phase velocity u/k. The problem encountered in the quantum case is 
that the numerators are not linear in k, so that their fc-dependences are not removed by 
this division. The problem encountered in the collisional case of [TUl [16] , where uj in the 
denominator of the ion term is replaced by u + z/j, where v x is the ion- neutral collision 
frequency, is that v\jk is not independent of k (unless v\ is assumed proportional to k) so 
that the denominator no longer depends only on u/k. Similarly, in the magnetised case 
the cyclotron frequencies, u cs , would complicate the fc-dependences of the denominators. 
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Figure 1. Coordinate system used for numerical evaluation of the solid-angle integral 
in @. 



3. Numerical formulation 

Two systems of coordinates for evaluating the solid-angle integral in Ml) numerically 
suggest themselves: spherical polars with axis along vo, and spherical polars with axis 
along r as shown in figure [TJ We adopt the second choice because it appears to have 
the advantage that it is the natural choice when deriving the bare potential term in (J6J), 
where the singularity at k-r = is of particular concern and v does not appear. In 
this coordinate system we have 

k-v = f/,v\\ + v/l — /U 2 cos0fj_ , (10) 



where v\\ = v -r = uocos^, v± = |v Xr| = i^sin^, and /j, = k-r. 

In these coordinates the element of solid angle in Q is given by dVt = dfxdcj), with 
the ranges of integration being \i e [— 1, 1], e [0, 2tt]. By using the reality condition, 
the ranges of the jx and 4> integrations can be reduced by half, so (J6j) becomes 



(p(r,ij)) 



2Mr) 



n 



71 

2 + 



dxlm 



$ T] (xV$>) 



with 



$ I w cos ip — h v o sin ip cos ( 
r 



1/2N 



1/2 



(11) 



(12) 



The first term in the integrand of the ^-integral (which arises from the delta-function 
in the Plemelj formula) gives the bare Coulomb potential ipo- The second term, the 
integral over x = rfi, gives the dressing potential from the screening cloud. As will be 
commented on further below, these two terms almost cancel in the far-field of Debye- 
screened regions. 
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Figure 2. Visualisations of the wake of a superthermal but subsonic test particle 
moving at Vg = 4,v t hi = 4w p j/£;Di i n a Maxwellian-ion plasma with kjy c = (colour 
online), (a) Contour plot of ip(r) / tp {r) . This figure agrees with figure 13 of [7J, the 
dot-dashed lines delineating the "thermal Mach cone." (b) VQ(p(r t n) vs. r/v behind 
the particle. This figure agrees with figure 12 of [7J. (c) tp(r, if)/<Po(r) vs. r/v behind 
the particle. This is a transect of (a) along the negative y = axis, extended to 
x/v Q = -20. 



4. Numerical results: validation against MSc plots 

We have reused the subroutine ETA listed in reference [5] (a listing with a README file 
and a test program being provided online as supplementary material for this paper) for 
calculating the special function rj(z), (A.l|), in a FORTRAN 77 reconstruction of the 



original FORTRAN program for computing cp(r,ip) from (11). This code computes the 
double integral over and x using Romberg integration [25] with the maximum size 
of the approximation matrix set to 10 x 10 in both integrations. The Free Software 
Foundation's gcc-based gf ortran [2TJ was used to compile and link the program. 

In the calculations we used units used such that fcoi = w p ; = 1 and q = 47r£ so that 
(po(r) = 1/r. Assuming the particle velocity to be much less than the mean electron 
velocity, we approximated the electron polarisation function <3> e with its static value, 
$ e (0) = k^, e and used the dawson routine [25J to calculate the ion polarisation function 
from g. 

Figures |2^a) and |3^a) show contour plots of (p(r)/(p (r) in the x/vo,y/vo plane, 
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Figure 3. Visualisations of the wake of a supersonic test particle moving at 
vq = 32uj pl /km = 4C S in a Maxwellian-ion plasma with ^Do = km/ 8. (a) Contour plot 
of (p(r)/ipo(r). (b) ip(r,ip)/ipo(r) vs. ip (in degrees) on the small arc r — 1.5«o/fcDi) on 
the right of (a). This figure agrees with figure 14 of [7J. (c) f/(flo on the middle arc, 
r = 7.5wo/fcDi- (d) ^/^o on the left-most arc, r = 15i>o/fcDi- 



where x and y are such that r = xx + yy + -zz, with the unit vector x in the direction 
of Vq. The other panels show the behaviour of <p along transects as described in the 
captions. The figures are both for test particles moving faster than the ion thermal 
speed, figure |2] showing a subsonic case and figure [3] showing a supersonic case (see 
captions) . 

These two figures were produced by interpolation from 100 x 100 grids above the 
x-axis, with the contours below the axis being obtained by reflection. The 10,000 
evaluations each of the Fourier inversions ^ were calculated as described above. Using 
gf ortran on a MacBook Pro, the scan for figure [2] took 311 seconds, while that for 
figure [3] took 942 seconds, giving an average of 0.03 s and 0.09 s per Fourier inversion, 
respectively. The interpolations and plots were done with Mathematica [26J. 

In figure [3^a) the disturbance behind the particle is seen to be limited to a Mach 
cone (dashed lines) of half width about 15 degrees. Note that the Mach cone does not 
contain a shock wave owing to the dispersive character of ion acoustic waves. Asymptotic 
analyses performed in reference [5] were published in [7] . Further analyses, including an 
interpretation of the wave structure observed in figure [3] based on ray tracing, will be 
published elsewhere. Here we limit ourselves to a qualitative heuristic interpretation of 
the wake. 

In the rest frame of the plasma, the screening response excited by the bare potential 
of a test particle as it passes a given point occurs with a time lag, producing a region 
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of opposite potential immediately behind the particle but travelling with it, as shown 
in figure [3](b). In figure |3^c,d) it is seen that the sign of this near-field disturbance is 
preserved away from the x axis in the potential minimum just inside the Mach cone in 
the far field. This is a wave train consisting of a superposition of ion acoustic waves with 
phase velocities given by kk-v , diminishing in amplitude away from the test particle 
due both to spreading and Landau damping. 

In the Debye-screened regions outside the Mach cone there is strong cancellation 



between the bare potential and screening cloud terms in (11). Thus, when ip is very 
small the relative error can become significant, leading to some numerical "noise" causing 
spurious zero contours in the upper and lower right-hand corners of figure [3j This was 
filtered out by setting <p/<po to 0.01 when it fell below 0.01 in absolute value outside 
the Mach cone, but it may be possible to eliminate the problem by using Q instead of 
(JgI) and deforming the contour of the numerical /i-integration away from origin to avoid 
the singularity there, thus avoiding the split between bare and screening potentials and 
handling Debye cancellations analytically. 

5. Quasiparticles 

As indicated in the Introduction, the original motivation of this work [5, Chapters 2 
and 3] was mainly to achieve a better understanding of plasma kinetic theory via a 
visualisation of dressed test particles. 

Another motivation was to test the calculation by Pines and Bohm [27] of dynamical 
screening using a classical version of their quantum-mechanical collective coordinate 
approach j28l[29]. Their calculation in fact failed this test JTJ p. 9], but more recently 
a corrected quantum-mechanical collective coordinate approach has been developed to 
calculate the interaction of dust particles in plasmas [30J . 

While Chapter 4 (unpublished) of [5] was on a quantum field theory treatment of 
electron-photon scattering in a plasma, including a test particle calculation, the first 
author has not used quantum mechanics per se since, except for a short project on 
statistical mechanics of a thin film [31] after submitting his MSc thesis and before 
commencing PhD studies at Princeton. 

However his experience with the quantum approach for calculating nonlinear wave- 
particle and wave-wave interactions led naturally to the thought that perhaps the 
relative ease with which this is done in quantum field theory [321 [331 e -g-] is due to 
the power of the Lagrangian and Hamiltonian formalisms developed in the field, rather 
than being due to anything intrinsic to quantum physics. Perhaps, if similar effort 
were applied to developing Lagrangian and Hamiltonian methods for classical plasmas, 
similarly powerful results would follow. This thought, combined with the one related 
to the dressed test particle picture described below, has been expressed in much of his 
later work. 

In quantum field theory, especially as applied in condensed matter physics, dressed 
test particles are regarded as "quasiparticles", related to bare particles but with 
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properties altered by the interactions between the bare particles. In quantum theory the 
first step is often to write down a Hamiltonian operator consisting of an unperturbed 
(bare) part and an interaction term, and then to diagonalise the Hamiltonian using a 
Bogolyubov transformation. The analogue of this in classical Hamiltonian theory would 
be to find a canonical transformation to a normal form in which explicit many-body 
interactions were reduced to an irreducible residual part, with the collective interactions 
being incorporated in a renormalised "unperturbed" Hamiltonian. 

The Pines-Bohm [27] development of collective coordinate theory for particles 
interacting via Coulomb potentials was very much based on classical canonical 
transformation theory. However, the first author began to develop a quasiparticle 
formalism not in a many-body Hamiltonian context, but in a Lagrangian approach 
to wave-background interaction in an ideal magnetohydrodynamic (MHD) fluid [34J, 
which showed that the conservation of wave action arises as naturally in classical 
mechanics as in quantum mechanics. However, another expression of the quasiparticle 
concept was in developing an oscillation-centre theory of turbulent phase-space diffusion 
[35] |36| I3~7] . which used canonical Hamiltonian perturbation theory. The operator 
method for classical canonical transformations introduced in [37] triggered the adoption 
of Lie methods as an important tool in theoretical plasma physics [38] . 

As it became apparent that the separation of the transformed Hamiltonian into a 
renormalised, non-interacting part and a residual interaction part was highly non-trivial 
and intimately related to developments in Kolmogorov-Arnold-Moser (KAM) theory, 
more recent development of this line of thought has been in papers aimed at developing 
a canonical theory of almost-invariant tori in 1 1 Hamiltonian systems via a pseudo-orbit 
approach [331 EH El E2] , with physical application to the description of magnetic fields 
in non-axisymmetric toroidal plasma confinement systems [431 HH 05} SHI H7J. Whether 
this approach will ever succeed in finding a classical Hamiltonian derivation of the 
Thompson-Hubbard dressed test particle picture is yet to be seen, but the development 
of modern Poisson-bracket perturbation formalisms [H] may point the way. Also, recent 
developments in the collective coordinate approach [3D] look very promising. 

Other concepts picked up during the first author's MSc research, such as the use 
of special functions, asymptotic expansions and numerical analysis also of course have 
found expression in various ways in his subsequent publications but we do not attempt 
to analyse these here. 

6. Conclusion 

We have presented analytical and numerical details of a novel approach, introduced in 
reference [5], to inverting the Fourier transform of the screened field of a test particle, 
and have illustrated the method with results from two cases, subsonic and supersonic 
test particles. 

The lack of computational detail in the literature, and the rapid increase in 
computer speed over the years, makes direct comparison with other approaches difficult. 
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However, the approach presented here, being based on classical analysis techniques for 
special functions and sound numerical analysis principles, is expected to be relatively 
fast. Anecdotal evidence appears to bear this out. The reason that the method has 
not been presented before is presumably that the "auxiliary function for sine and cosine 
integrals," f(z) [19], which could have been used instead of our function rj(z), is not 
well known in the physics community, and no subroutines for its rapid evaluation have 
previously been available. A primary purpose of this paper is to increase knowledge of 
this function and perhaps stimulate the development of better subroutines than the one 
we have provided online (see the online README file for error plots). 

The use of this special function imposes certain restrictions on the applicability of 
the method, making it inapplicable to quantum, magnetised and collisional plasmas as 
discussed in section 12. 1[ 

In section [5] we have also sketched how the dressed-test-particle picture is related 
to the concepts of quasiparticles, oscillation centres and pseudo-orbits. This section has 
been included to help relate the present paper to others in the cluster of papers for 
which it is intended. 
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Figure 4. (a) Real and (b) imaginary parts of rj{x + iy) showing their respective 
antisymmetry/symmetry ( |A.9[ ) under reflection in the imaginary axis of the (x + iy)- 
plane, cut along the negative imaginary axis. 



Appendix A. The function r] 



The function f](z), first introduced in reference [3], was used in both the analytical 
and numerical work in the unpublished thesis. Much of reference [5] was published 
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in reference [7], but Appendix I, giving details of the mathematical properties and 
numerical calculation of i](z), was omitted. Thus we reproduce the contents of 
Appendix I below, including a small amount of additional material for clarity. 

The function is closely related to the exponential integral [49, §5.1.1], and may be 
computed from tables of this function. As shown below it is even more closely related to 
the auxiliary function for sine and cosine integrals, f(z) [121 §5.2.6] and [50], but does 
not appear to have been defined before so the notation is our own. 

Define rj(z) by 

/ \ z r° , e " . vr , . , . 

V{z) = ~ / dx ^~, — a for \ ar Sz\<-, (A.l) 

% Jo ar + z l 2 

and analytically continue the function into the left half plane, cutting the complex plane 

along the negative imaginary axis (see figure El). 

It may be shown that r\{z) has the alternative integral representation 

™ e~ l 
h z 2 — t 2 
The form most useful for our purposes is the identity 
r°° expiax 
i Jo x 2 + (3 2 

For instance, we can relate ■!](■) to the auxiliary function f(z) [HI §5.2.12] by putting 
(3 = 1, a = iz to give f(z) = irj(iz) for Re z > 0. Then, as f(z) is defined on the 
complex z-plane cut along the negative real axis [121 §§5.2.2, 5.2.6], replacing z with 
— iz we have, by analytic continuation, an alternative definition for rj(z), 

-77" i7T 

rj{z) = -if(-iz) for - - < Arg^ < — , (A.4) 



7](z) = z J dt — for < argz < 7r . (A. 2) 



poo p"Y"T) %OlOC 

/ o no dx = rj(a/3) for Imo > 0, Re/3 > . (A.3) 

Jo x 2 + V 



the cut now being along the negative imaginary axis as stated after (A.l) 



As will be found useful in the subroutine to be described below, we may also show 



from \k.2\ that 

1 



<+)f~\ „-*p(-)/ 



V (z) = -[e*E?>(z)-e-*Er>(-z)], (A.5) 

where E{ (z) are analytic continuations of the exponential integral E\(z) [49, §5.1.1] 
and [51] to the z-planes cut along the negative /positive imaginary axis (analytic in the 



upper/lower half planes), respectively [see below (A. 6)]. 
These functions can be represented as, 



-zy 

nn\ 



E{ ± \z)= - 7 -l n C±)(z)-£ 

n=l 

= -7-ln (±) (^) + Ein(z) (A.6) 

where 7 is Euler's constant, Ein(z) [51] is an entire function, and In' \z) = ln\z\ + 
i[arg(^f iz) ± tc/2] = \nz + i[arg(^f iz) — argz ± tt/2]. Thus, comparing ( |A.6 ) with 
reference [491 §5.1.1], 

E ( { k: \z) = E x (z) - i[axg(Tiz) - arg(z) ± tt/2] . (A.7) 
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As hr + '(z) and hr~'(—z) are both denned on the complex plane cut along the 
negative imaginary axis they can differ only by a constant: hr -*(— z) = ln^ + \z) — in. 
In (A. 5) this gives 

rj(z)= -[-f + \n {+ \z)]smhz-—e~ z + 1 [e z Em(z)-e- z Em(-z)} (A.8) 

= - 7 + ln w (z) smhz e~ 2 + 2 H 1 h U 7 , 

L/ v ;J 2 36 7200 V y ' 

from which can be seen that T](z) still has a branch point at the origin, but, unlike E\(z), 

it remains finite there. 

The following symmetry properties, reflection about the imaginary axis and 

reflection about the real axis (with exponential correction), can be proved from (A.8) 

and may be used to continue rj(z) out of any quadrant into the other three quadrants, 

r,(z)= -v(-z*)\ (A.9) 

v (z) = r)(z*)* - me~ zsvez , (A.10) 

where 

sre z = sgn (Re z) (A. 11) 

and * denotes complex conjugation. These symmetries are apparent in the plots in 
figure |4j 

The asymptotic expansion for large \z\ may be obtained from that for f(z) |^9l 
§5.2.34] 

1 / 2' 4' \ 7T 37T 

r/(z) ~ - ll + — + — + ... J as |z|^oo, -- < Arg2; < — .(A.12 

However, the successive terms in this expansion do not decrease sufficiently rapidly for 
it to be very useful for computational purposes over the range of \z\ in which we are 
interested. 

All the derivatives r]( n '(z) of rj{z) may be expressed in terms of rj{z) and r)'(z), 

,W(») =»W-^-^--..-™; «-en, (A.13) 

,W( 2 ) = V( 2 ) + ^^ + ^^ + ... + ^; nodd. (A.14) 

The two relations above are used to calculate r](z) for intermediate \z\ by extrapolating 
from known values of r](z) and rf(z) using Simpson's rule. The known values are 
calculated from table 5.6 of reference (49J using (A. 5) and (A. 10) and the relations: 

V , (z) = 1 -[e z E[ + \z) + e~ z E[-\z)] (A.15) 

rj'(z) = n'(z*)* + sre(z)me- ZSICZ , (A.16) 



where sre(-) is defined in (A. 11) and 

E[ ± \±z) = E 1 (±z) for lmz>0, (A. 17) 

E^zf = E 1 (z*) . (A.18) 
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For large \z\, Laguerre integration [321 §25.4.45] of (A. 2) is used, 

z A — xt 



r] {z) = zJ2 2 _' 2 +Rn] ImoO, (A.19) 



=i 

\2 



Kl < TTZ^ZT , (A-20) 



(Imz] 



2n+l 



with the abscissas Xi and weight factors W{ obtained from |1H1 table 25.9]. The error 
bound is pessimistic at large |Rez|, as this equation provides accuracy of six significant 
figures, even for z real, if \Kez\ > 16. 

The FORTRAN IV subroutine developed for efficient calculation of rj(z) using the 
methods described in this appendix was listed in reference [5j. It is accurate to around 2 
parts in 10 4 and evaluates r)(z) in about 0.5/is, so this innermost integral over < k < oo 
in the Fourier inversion may be regarded as performed analytically, leaving only the 2- 
dimensional integral over solid angle to be done numerically. 

To produce the numerical results presented above we used this core subroutine in 
our reconstruction of the original screened field program (a testament to the backwards 
compatibility of FORTRAN). For the historical record we provide the FORTRAN source 
code, and a README file describing its use, as supplementary data in the online version 
of this paper. However, for accurate calculations it would probably be preferable to use 
the professionally written FORTRAN code for E\(z) developed by Amos [52"] . 
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